Apparent phase transitions in finite one-dimensional sine-Gordon lattices 
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Illes Balears 



We study the one-dimensional sine-Gordon model as a prototype of roughening phenomena. In 
spite of the fact that it has been recently proven that this model can not have any phase transition [J. 
A. Cuesta and A. Sanchez, J. Phys. A 35, 2373 (2002)], Langevin as well as Monte Carlo simulations 
strongly suggest the existence of a finite temperature separating a flat from a rough phase. We 
explain this result by means of the transfer operator formalism and show as a consequence that 
sine-Gordon lattices of any practically achievable size will exhibit this apparent phase transition at 
unexpectedly large temperatures. 

PACS numbers: 05.70.Fh, 05.10.-a, 68.35. Ct, 64.60. Cn 
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I. INTRODUCTION 

More than fifty years ago, van Hove |j| proved that 
true thermodynamic phase transitions, defined as singu- 
larities of the free energy, could not occur in a class of 
one-dimensional (ID) systems, a result later extended to 
lattice systems in the same class by Ruelle In spite 
of the fact that the conditions for van Hove's theorem to 
apply were clearly stated Q (see also Q), there is nowa- 
days a very general belief that ID systems cannot exhibit 
phase transitions unless they have long range interac- 
tions. This misinterpretation of van Hove's mathemati- 
cal results has been reinforced by the abuse of Landau's 
argument about the entropic contribution of domain 
walls to the free energy. This argument, being physically 
very intuitive and useful, is not a rigorous result (assump- 
tions such as a dilute concentration of domain walls are 
made along the way) and, furthermore, it does not apply 
to every ID system. In fact, there are many examples of 
ID systems with true thermodynamic phase transitions 
H 0- @- 1[ ||] wn i cn i unfortunately, have remained largely 
unnoticed. 

In the more specific context of models of growth pro- 
cesses 0, pLltl , the unsustained belief discussed above 
is often translated by saying that ID interfaces are al- 
ways rough. This is actually not the case, as shown in 
the early eighties with several examples J?], || . Only re- 
cently, two of us JT^ ] have proven a theorem showing that 
a wide family of ID models, including the sine-Gordon 
model as a particular example, cannot have phase tran- 
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sitions, this being, to our knowledge, the first time that 
a rigorous proof of that kind has been given (some non- 
rigorous, phenomenological arguments in the same direc- 
tion had been proposed earlier ]T^|). However, as we will 
see below, this theorem turns out to be in conflict with 
some numerical simulation results which seem to provide 
strong evidence supporting the existence of a roughening 
phase transition in the ID sine-Gordon model. In view 
of the fact that simulations are very often the only way 
of studying a large class of models, it is most important 
to understand this contradiction in order to distinguish 
between true and apparent phase transitions. 



To the above end, in this paper we focus on the ID 
sine- Gordon model as a canonical example, widely appli- 
cable and representative of the phenomenology of many 
model systems juj (see also |lf| for a review). Thus, 
in Section [0] we give results of simulations that suggest 
the existence of a phase transition at a (not necessar- 
ily small) nonzero temperature. By means of a transfer 
operator approach and using the probabilistic meaning 
of the corresponding eigenfunctions 



III 



we 



in Section 

analyze the origin of this behavior; from this analysis 
we are able to conclude that such apparent phase transi- 
tions will occur not only for lattice sizes achievable within 
the present computational capabilities, but also for very 
much larger lattices. Finally, in Section IV we discuss the 



consequences of this result, which we believe are relevant 
for computational studies where no analytical support 
exists. Furthermore, additional important implications 
of our research for experimental studies of small systems 
far from the thermodynamic limit are also considered. 
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II. NUMERICAL SIMULATIONS 

The ID sine-Gordon model is denned by the following 
Hamiltonian: 



H 



N T 



h,) 2 + V [l-cos(h l )}}, (1) 



where N is the number of lattice nodes (or the system 
size), J is the coupling constant and — oo < hi < oo is 
a real variable on site i. We assume periodic boundary 
conditions ho = /ijv- For visualization of our results, we 
interpret hi as the height of a surface above site i of the 
lattice; then, the two terms of the Hamiltonian corre- 
spond, respectively, to surface tension and to a local po- 
tential (of strength Vq) favoring multiple values of 2ir for 
the height, representing that growth takes place prefer- 
entially by addition of discrete units (layers). For surface 
growth on two-dimensional (2D) substrate lattices, this 
interpretation has proven itself rather fruitful in the past 
(see [l7], [l9| and references therein). However, we 
want to stress that the results we present in this paper 
are independent of any specific interpretation one makes 
of hi. In fact, previous studies of the ID sine-Gordon lat- 
tice jl4| 20 were more interested in understanding the 
role of solitons in statistical mechanics models |l5| than 
in any particular application. 

In our study, we concentrate on two magnitudes in 
order to characterize the model behavior: the surface 
width or roughness, 



10 



where 



1 N 

-Yh 



(2) 



(3) 



is the mean height, and the height-difference correlation 
function, 



C(r) 



' 1 N 



(4) 



Averages (■ • • ) are to be understood with respect to a 
statistical weight given by the Gibbs factor, e~ u / T , at 
equilibrium at a temperature T. 

The above defined are crucial quantities in the 2D ver- 
sion of the model. This exhibits a Kosterlitz-Thouless 
type phase transition from a low temperature flat phase 
to a high temperature rough phase |l7], ^l], ^2], p3| . 
In the flat phase, small systems have a size dependent 
width whereas the width of large systems is independent 
of the size. The crossover system size separating both 
regimes is closely related to the correlation length, which 
is finite in the low temperature phase, and can be defined 



as the distance beyond which the height-difference corre- 
lation function saturates. On the contrary, in the rough 
phase the correlation length is infinite and, correspond- 
ingly, the roughness increases (logarithmically in the 2D 
case) with the system size for all sizes, i.e., it is also in- 
finite in the thermodynamic limit. In ID, the theorem 
proved in [ fl2| prohibits any phase transition, and at all 
nonzero temperatures the system is in the rough phase, 
the roughness increasing linearly with the system size. 

In the lack of detailed analytical results, the statistical 
averages can be computed approximately by means of nu- 
merical simulations. This kind of analysis has become a 
routine tool for the study of the equilibrium properties of 
many models and a problem of interest is to extract from 
the numerical studies, necessarily performed in finite size 
lattices, the asymptotic behavior in the thermodynamic 
limit. We will show that a naive extrapolation of the fi- 
nite size results for the ID sine-Gordon model can lead 
to erroneous results concerning the existence of a phase 
transition at a finite value of the temperature. 

For our numerical study, and in order to assess the va- 
lidity of our results, we have used two completely differ- 
ent procedures: Langevin dynamics and parallel temper- 
ing Monte Carlo. The Langevin dynamics procedure has 
been widely used in this context with very good results 
]Tt| , [Ts| , pj| , and it consists in the numerical integration 
of the Langevin equation following from the Hamiltonian 
H: 



dh t {t) 
dt 



on 

dhi(t) 



■Tk(t), 



(5) 



where rji are Gaussian white noises obeying the 
fluctuation-dissipation theorem at temperature T, i.e., 



(vj(t')m(t)) = 7r8ijS{t--e). 



(6) 



A major problem with the Langevin dynamics is the pres- 
ence of systematic errors, in addition to the unavoidable 
statistical errors, due to the finiteness of the time step 
used in the numerical integration (in our studies, we have 
used a stochastic Heun method |24|]). 

The second procedure, parallel tempering, is the one 
we have mostly relied on. The reason is that it is a 
very efficient algorithm to prevent the system from being 
trapped in local minimum energy configurations. Par- 
allel tempering requires any Monte Carlo method which 
generates representative configurations at a given tem- 
perature. In this case, we have implemented a heat bath 
algorithm p5[ , in which new values h\ for the height at 
site i are proposed according to the rule 



h' = 



ha 



(7) 



£ being a Gaussian random variable of zero mean and 
unit variance, and are accepted with a probability 
min[l, e- SH / T ] with SH = -V [cos(^) - cos(/i*)]. The 



parallel tempering algorithm considers then simultane- 
ous copies of the system at different temperatures, al- 
lowing exchange of configurations between them. The 
exchange occurs after enough configurations have been 
generated at each temperature for a time greater than 
the energy autocorrelation time (see, e.g., [^6[ ^7) for 
details). This is particularly efficient for low tempera- 
ture configurations, which are most susceptible to being 
trapped in metastable regions. All the results reported 
in this paper have been obtained by means of this paral- 
lel tempering Monte Carlo algorithm, although we have 
checked that Langevin dynamics produces the same re- 
sults (quantitatively within error bars). 

Before proceeding with the discussion of the simulation 
results, let us recall the theoretical background. As we 
have already mentioned, the theorem proven in |12[] im- 
plies that the system is in the rough phase at all nonzero 
temperatures. This can be interpreted in terms of general 
renormalization group arguments (for a renormalization 
group study of the 2D sine-Gordon model, see ]2l], |||, p3| ; 
see also jnj, [ll| for a more general perspective): in the 
ID model, fluctuations should be enough to effectively 
suppress the potential part of the Hamiltonian at any 
nonzero temperature, leaving as the only relevant term 
the discrete gradient (surface tension). In that case, the 
sine-Gordon model behaves effectively as the Edwards- 
Wilkinson model p8| , whose associated Langevin equa- 
tion (|^) is simply the discrete linear diffusion equation 
with additive noise, and all the properties of interest can 
be calculated. This approach has been very successful in 
characterizing the 2D sine-Gordon model behavior Jl9| , 
and particularly in locating the roughening transition 
temperature. In our ID case, it is easy to show that in the 
Edwards- Wilkinson regime the roughness must scale lin- 
early with the system size for nonzero temperature ]29[ ] . 

Figure [I] displays the simulation results for the rough- 
ness for several system sizes as a function of temperature. 
For the sake of definiteness we have chosen J = Vq = 1 
in Eq. ([!]); other choices yield the same qualitative re- 
sults. Figure |l| shows that, as expected, the data tend 
asymptotically, for high values of the temperature, to 
the Edwards- Wilkinson result, u 2 /N = T/12 (29). This 
linear scaling of the roughness with system size indi- 
cates clearly that the surface is rough at high temper- 
atures. However, the main plot in the figure indicates a 
clear change of behavior around a temperature T ~ 1. 
In fact, as shown in the inset zoom, at low tempera- 
tures (T < 0.8) the linear scaling dependence of the 
roughness with system size is lost and in that region the 
roughness becomes fairly independent of the system size, 
a behavior that according to the preceding discussion 
would correspond to a flat phase. To obtain a theoret- 
ical prediction for low temperatures, we have analyzed 
yet another linear model, in which the cosine term in 
the Hamiltonian is substituted by a parabolic potential, 
Vb(l — cos(hij) — ► Vo/if/2. Such a model is flat at all 
temperatures (basically because the parabolic potential 
confines the surface to lie around its minimum) and, as 
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FIG. 1: Roughness scaled by the system size, u 2 /N, vs tem- 
perature for the ID sine-Gordon model. System sizes are as 
indicated in the figure. The dotted line corresponds to the 
theoretical value uj 2 /N — T/12 obtained for the Edwards- 
Wilkinson model. Inset: Zoom of the low temperature region 
showing the lack of scaling. The dashed line corresponds to 
the parabolical approximation discussed in the text, and is 
given by uj 2 = T/^/V 2 +4V . 
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FIG. 2: Log-log plot of the height-difference correlation func- 
tion scaled by temperature vs distance scaled by the system 
size. Bottom to top, temperatures are 0.0956, 0.2407, 0.7029, 
0.8115, 0.9896, 1.1689, 1.4819, 1.9016, 2.5562, 3.7044, 6. Also 
plotted are the predictions of the Edwards- Wilkinson model 
and of a parabolic approximation (see text). Error bars are 
typically as shown in one of the curves. 



can be observed in Fig. [I], agrees very well with the nu- 
merical simulations at low temperatures. We thus see 
that, in spite of the fact that we know that no phase 
transition can take place in this model, something very 
similar to a phase transition from a flat to a rough phase 
appears in the simulations for a temperature T* ~ 0.8. 

The above indications arc reinforced by looking at the 
correlation function, shown in Fig. ||, where we can see 
that for low temperatures there is a finite correlation 
length, whereas for high temperatures the correlation ex- 
tends as far as half the system size (recall that we use pe- 
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riodic boundary conditions). It can also be appreciated 
again that for high temperatures the simulation results 
reproduce quite well the Edwards- Wilkinson prediction. 
In the opposite limit, much as it occurs with the rough- 
ness, the correlation behaves like the parabolic approx- 
imation. We have also studied other magnitudes, such 
as the specific heat, finding exactly the same result as 
Schneider and Stoll |Q that the specific heat exhibits a 
maximum at a value somewhat higher than T* . We want 
to stress that all this is very reminiscent of the behavior 
in 2D, where the maximum of the specific heat is inter- 
preted as a Schottky anomaly (see (ll| and references 
therein). Hence, from the available numerical evidence, 
we would be forced to conclude that there is a roughen- 
ing transition at a temperature T* in the ID sine-Gordon 
model were not for the theorem in fL2l. 



III. TRANSFER OPERATOR APPROACH 

In order to understand the numerical results, we will 
use the transfer operator approach jl], ||, |l4|, specif- 
ically, to make the connection with the presentation in 
fl2f , we rewrite our Hamiltonian ([!]) rescaling the hi vari- 
ables by a factor 2ir, i.e., 



N 



.An 2 J 



H = {-o-fai-i - h t f + V [l - cob(2tt/h)]}. (8) 

i=i 

In this notation, the preferred values for hi are the in- 
teger numbers, and the Hamiltonian is invariant under 
the transformation hi i— > hi + 1. This means that we 
can choose for convenience hi £ (—5, 5] without loss of 
generality. With this choice in mind, the corresponding 
partition function becomes 

/1/2 poo />oo 

dhi / dh 2 --- / dh N e- 0H , (9) 
-1/2 J-oo J -co 

(3 being the inverse temperature in units of the Boltz- 
mann constant. We now write hi = 7ij + 4>i, with m G Z 
and — i < <pi < I , with i = 1, . . . , JV. Let us define 

00 

V(fi,<f>,0)= e-^^^+^V 1 ™ 9 , (10) 



n— — 00 



a 27r-periodic function of 9, and the operator 
,1/2 

T ,ef(<t>) = / ' 
J -1/2 



b'T , e (4>,^)f (<!>'), (11) 



Tp fi {<l>,tf) = V{f3,4>-ct>',e) x 

0, 



(12) 



x exp j-^V [2 - cos(2tt0) - cos(2tt(//)] j ; 

with these definitions, the partition function can be writ- 
ten as 

1 r 



27r J-7T 



(13) 



T p t e is called the transfer operator for this model. Using 
the operator properties, it can be shown (see (I2) for 
details) that in the thermodynamic limit (JV — > 00) 

TrfT^r = 5>„(/?,tf)r 

n>l 

= m(P,9)[X max (P, 6)] N [l + o(l)} (14) 

where A„ are the operator eigenvalues, necessarily real 
and isolated, A max is the maximum eigenvalue, necessar- 
ily positive, and m(/3, 9) is its multiplicity, necessarily 
finite. Finally, Laplace's method yields the free energy 
in the thermodynamic limit: 

-/?/(/?)= lim ±lnZ N ({J)= max ln[A max (/3, 9)). 

JV— *00 iV — 7T<y<7T 

_ (15) 
Based on this expression, in |12] it was proven that 
the maximum of A max (/3, 9) occurs at 9 = 0, and that 
A m ax(/3, 0) is analytic for /? > 0, which leads to the con- 
clusion that the free energy itself is analytic for all j3 > 
and, subsequently, that there are no phase transitions in 
the model. However, that is not all the information we 
can obtain from this approach, as we will now show. 

In |1^| , the squared modulus of the eigenfunction of 
the largest eigenvalue is interpreted as the probability 
density for the hi variables. As the transfer operator in 
JTij is different from the one we are using here, and for 
the sake of completeness, we now proceed to show that 
we can resort to the same interpretation here. Leaving 
out irrelevant constants, and keeping in mind that as we 
have just said the only contribution to the free energy 
comes from 9 = 0, we can write 



N 



Z N (f3) = TrfT^o) 

We can now compute averages of functions g(4>j) 
following way O] : 



(ff(fc)> 



N-j 



TrT 



18,0 



N 



(16) 
in the 



(17) 



In the thermodynamic limit, it can be shown that the 
previous equation becomes 



1/2 

(gib)) = I d</>M</>)| 2 <?W>) 
1/2 



(18) 



where tpo((f)) is the eigenfunction of the largest eigenvalue. 
We thus see that, indeed, |y>o(0)| 2 can be understood as 
the probability density for the height to be in the interval 
(—5,5] modulo 1. 

In order to apply this result, we must compute the 
eigenvalues and eigenfunctions of the transfer operator. 
This has to be done numerically: to this end, one has to 
discretize the operator and transform it into a matrix (see 
|j, 14 , sec 1 30 for a detailed account). The advantage of 
the present formulation of the transfer formalism is that 
integrals are to be carried out on a finite interval, and 
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FIG. 3: Kink density as estimated from the probability that 
the value of the height goes over the potential maxima at <f> = 
±1/2 (see text). Horizontal lines correspond to the inverses 
of the system sizes studied. For comparison, the approximate 
result obtained by Schneider and Stoll jl3 is included as a 
dotted line. Inset: squared modulus of the eigenfunction of 
the largest eigenvalue at low temperature. 
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FIG. 4: Roughness vs temperature for different system sizes 
(right to left: N = 50, 125, 250, 500, 1000, 2000). The arrows 
mark the temperature predicted by our criterion based on the 
eigenfunction for each of those system sizes. Inset: estimates 
for T* vs 1/lniV. Points are obtained from our criterion, line 
is a linear regression fit. 



therefore we do not need to introduce ad hoc any cutoff 
as in the case of integrals on an infinite interval, thus 
eliminating one possible source of error or inaccuracy. 
Specifically, in our numerical diagonalization procedure 
we have used 2001x2001 matrices; we have checked that 
increasing their size does not change the results signifi- 
cantly. As another check, we have computed the specific 
heat from the numerically computed eigenvalues, finding 
perfect agreement with the output of the simulations. 

The inset in Fig. || shows the squared modulus of a 
typical eigenfunction at low temperature. The interpre- 
tation in terms of probability density indicates that prob- 
able values for the height lie close to the minimum of 
the potential, i.e., to = 0, whereas values close to the 
maxima of the potential at (f> = ±1/2 are very unlikely 
to occur. We can associate the probability of taking a 
value of <fi — ±1/2 with the probability of formation of a 
kink or step, as once the height is at a maximum it can 
cross over to the neighboring potential well, thus giving 
rise to a kink. This interpretation suggests us to com- 
pare that estimate with the inverse of the system sizes 
studied in this work, which is a reasonable estimation 
of the probability of observing a kink in our numerical 
system. Figure [}] compares both quantities, making clear 
that for every system size there is a temperature at which 
the probability of formation of one kink becomes smaller 
than the inverse of the system size. In fact, the proba- 
bility of formation of one kink decays extremely rapidly 
below (orders of magnitude, note the logarithmic scale) 
the crossing temperature. 

Following the discussion above, it is very natural then 
to associate that crossing temperature with T*, the tem- 
perature at which we observe the apparent phase transi- 
tion in our simulations. This is very well confirmed by 
Fig. 0, in which the roughness obtained from our numer- 



ical simulations is plotted along with the temperatures 
predicted by our criterion above. We note in passing 
that this result allows us to understand the reason why 
the crude parabolic approximation coincides so well with 
the numerical results at low temperatures: the hypoth- 
esis underlying the approximation is that kinks do not 
form and the whole surface lies close to a single potential 
minimum. Indeed, we have seen that the probability of 
kink formation at low temperatures is certainly negligi- 
ble. 

Summarizing, we have support for the criterion, as con- 
firmed by the comparison of the transfer eigenfunction 
approach to the results of the numerical simulations, that 
the apparent transition temperature T* coincides with 
the one that yields a kink formation probability smaller 
than the inverse of the system size. This criterion allows 
us to make the following quantitative prediction: for a 
system of size N, there is an apparent flat-rough phase 
transition at a temperature T* ~ 1/ In N, as shown in the 
inset of Fig. ||. This can be easily understood if we realize 
that as kink formation is an activated phenomenon, the 
kink density follows an Arrhenius type law (with correc- 
tions, see |Ti|]). The dependence of T* on N follows from 
our criterion by imposing the proportionality of the kink 
density and 1/N. We can now estimate the size of the 
system needed in order to observe the rough phase all the 
way down to any given temperature, T*. Taking as an 
example T* — 0.1, which is certainly not small, we find 
that lattices of the order of 10 30 sites are required to en- 
sure a reasonable chance that kinks are formed during the 
simulation and the rough phase is observed for T > T*; 
we would still find such an exceedingly large system in 
the flat phase for T < T*. We have thus shown clearly 
that systems of any practically achievable size will always 
exhibit an apparent phase transition at a temperature far 
from T = 0. 
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IV. CONCLUSIONS 

In this paper, we have studied analytically and nu- 
merically, by Langevin dynamics and mostly by parallel 
tempering Monte Carlo simulations, the ID sine-Gordon 
model. We have found in the simulations that there exists 
a temperature at which an apparent roughening phase 
transition takes place. We have shown that it is possi- 
ble to understand the contradiction of such phenomenon 
with the theorem that prohibits phase transitions in this 
model jl2| through the analysis of the eigenfunctions of 
the corresponding transfer operator. The interpretation 
of these functions as probability densities makes clear 
that lattices of any finite size will always show an appar- 
ent phase transition, because the probability that kinks 
are formed becomes negligible below certain tempera- 
ture. We have also seen that even in extremely large 
lattices the apparent transition occurs at temperatures 
far from zero and, in fact, T* ~ 1/lniV. 

The results summarized above are relevant in a much 
broader context, basically in two directions. First, our 
conclusion should be kept in mind when analyzing the 
outcome of numerical simulations of models about which 
there is little or none analytical information. Were not 
for the fact that we know that such a phase transi- 
tion is not possible, we would have concluded from our 
simulations that the ID sine-Gordon model presents a 
roughening phase transition. In fact, simulations for 
the 2D sine-Gordon model yield results very similar to 
those presented here , although in that case we have 
a true phase transition according to several approxi- 
mate calculations including renormalization group results 
fl6[ p2| , |2l| . It is important to realize that finite size 
analysis, which in principle could signal that the transi- 
tion goes to T — with system size, becomes question- 
able if the values for the apparent transition at the sizes 
amenable within computational capabilities are still very 



far from T = 0. In addition, in our model, approximate 
analytical results in the low temperature limit support 
the existence of the non-existent phase transition. There- 
fore, we conclude that one has to be extremely careful 
with claims of this kind. 

The second direction that our work points to is related 
to the very nature of phase transitions. True thermo- 
dynamic phase transitions, understood as singularities of 
the free energy or its derivatives, can only take place in 
the thermodynamic limit, and no such transitions occur 
in finite size systems. Hence, the apparent phase transi- 
tion we see in our simulation can indeed be thought of as 
a true transition in the context of finite systems. What is 
more important, similar phenomenology is bound to arise 
in small, mesoscopic systems, certainly far from any ther- 
modynamic limit one can think of. As systems of that 
scale become more and more relevant both for theoret- 
ical and for applied reasons, the question of the defini- 
tion and nature of phase transitions gains importance. 
In this respect, this work hints that non-thermodynamic 
transitions may well be physically existent, or, alterna- 
tively, that computations and results in the thermody- 
namic limit do not represent well the fate of large but 
finite systems, even of very large, mesoscopic ones. 
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